//
// Created by 王岩松 on 2020/3/8.
//

#include "matrix.h"

using namespace std;
matrix::matrix() {
    n = m = 0;
}
matrix::matrix(int _n, int _m) {
    n = _n;
    m = _m;
    A.resize(_n, vector<double>(_m,0));
}
matrix::matrix(int _n, int _m, double val) {
    n = _n;
    m = _m;
    A.resize(_n, vector<double>(_m,val));
}
matrix::matrix(vector<vector<double>> _A) {
    n = _A.size();
    m = _A[0].size();
    for(int i = 0; i < n; ++ i) {
        assert(_A[i].size() == m);
    }
    A = _A;
}
matrix matrix::appendRight(matrix a) {
    assert(n == a.n);
    matrix ans(n,m+a.m);
    for(int i = 0; i < n; ++ i) {
        for(int j = 0; j < m; ++ j) {
            ans.A[i][j] = A[i][j];
        }
        for(int j = 0; j < a.m; ++ j) {
            ans.A[i][j+m] = a.A[i][j];
        }
    }
    return ans;
}
matrix matrix::appendDown(matrix a) {
    assert(m == a.m);
    matrix ans(n+a.n,m);
    for(int j = 0; j < m; ++ j) {
        for(int i = 0; i < n; ++ i) {
            ans[i][j] = A[i][j];
        }
        for(int i = 0; i < a.n; ++ i) {
            ans[i+n][j] = a.A[i][j];
        }
    }
    return ans;
}
matrix matrix::transpose() {
    matrix ans(m,n);
    for(int i = 0; i < n; ++ i) {
        for(int j = 0; j < m; ++ j) {
            ans[j][i] = A[i][j];
        }
    }
    return ans;
}
matrix matrix::gauss() {
//    auto output = [&](vector<vector<double>> C) {
//        for(int i = 0; i < C.size(); ++ i) {
//            for(int j = 0; j < C[i].size(); ++ j) {
//                cout << C[i][j] << " ";
//            } cout << "\n";
//        }
//    };
    vector<vector<double>> B = A;
    int MN = min(m,n);
    for(int i = 0; i < MN; ++ i) {
        for(int j = i; j < n; ++ j) {
            if(sign(B[i][j]) != 0) {
                for(int k = 0; k < m; ++ k) {
                    swap(B[i][k],B[j][k]);
                }
                break;
            }
        }

        if(sign(B[i][i])== 0) continue;
        double tmp1 = B[i][i];
        for(int j = i; j < m; ++ j) {
            B[i][j] /= tmp1;
        }

        for(int j = 0; j < n; ++ j) if(j != i) {
                double tmp2 = B[j][i];
                for(int k = i; k < m; ++ k) {
                    B[j][k] -= tmp2 * B[i][k];
                    if( sign(B[j][k]) < 1e-15 ) B[j][k] = 0;
                }
            }
    }
    return matrix(B);
}
matrix matrix::inverse() {
    assert(n == m);
    vector<vector<double>> B(n, vector<double>(m,0));
    for(int i = 0; i < n; ++ i) B[i][i] = 1;
    matrix ans = (*this).appendRight(matrix(B)).gauss();;
    for(int i = 0; i < n; ++ i) {
        ans.A[i] = vector<double>(ans.A[i].end()-n,ans.A[i].end());
    }
    ans.n = ans.m = n;

    // free
    for(int i = 0; i < B.size(); ++ i) {
        std::vector<double>().swap(B[i]);
    }
    std::vector<std::vector<double>>().swap(B);

    return ans;
}
void matrix::Free() {
    for(int i = 0; i < A.size(); ++ i) {
        std::vector<double>().swap(A[i]);
    }
    std::vector<std::vector<double>>().swap(A);
}

//std::vector<double> operator + (std::vector<double> a, std::vector<double> b) {
//    assert(a.size() == b.size());
//    std::vector<double> res(a.size());
//    for(int i = 0; i < a.size(); ++ i)
//        res[i] = a[i] + b[i];
//    return res;
//}
//std::vector<double> operator + (std::vector<double> a, double b) {
//    std::vector<double> res(a.size());
//    for(int i = 0; i < a.size(); ++ i)
//        res[i] = a[i] + b;
//    return res;
//}
//std::vector<double> operator - (std::vector<double> a, std::vector<double> b) {
//    assert(a.size() == b.size());
//    std::vector<double> res(a.size());
//    for(int i = 0; i < a.size(); ++ i)
//        res[i] = a[i] - b[i];
//    return res;
//}
//std::vector<double> operator - (std::vector<double> a, double b) {
//    std::vector<double> res(a.size());
//    for(int i = 0; i < a.size(); ++ i)
//        res[i] = a[i] - b;
//    return res;
//}
//std::vector<double> operator * (std::vector<double> a, std::vector<double> b) {
//    assert(a.size() == b.size());
//    std::vector<double> res(a.size());
//    for(int i = 0; i < a.size(); ++ i)
//        res[i] = a[i] * b[i];
//    return res;
//}
//std::vector<double> operator * (std::vector<double> a, double b) {
//    std::vector<double> res(a.size());
//    for(int i = 0; i < a.size(); ++ i)
//        res[i] = a[i] * b;
//    return res;
//}
//std::vector<double> operator / (std::vector<double> a, std::vector<double> b) {
//    assert(a.size() == b.size());
//    std::vector<double> res(a.size());
//    for(int i = 0; i < a.size(); ++ i)
//        res[i] = a[i] / b[i];
//    return res;
//}
//std::vector<double> operator / (std::vector<double> a, double b) {
//    std::vector<double> res(a.size());
//    for(int i = 0; i < a.size(); ++ i)
//        res[i] = a[i] / b;
//    return res;
//}
